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Abstract 


The  report  presents  an  approach  for  a  general  laminated 
shell  geometry  describable  by  orthogonal  curvilinear 
coordinates .  The  theory  includes  a  through  the  thickness 
parabolic  distribution  of  transverse  shear  stress. 
Additionally,  a  simplified  approach  that  allows  large 
displacements  and  rotations  is  incorporated.  The  theory  is 
cast  into  a  displacement  based  finite  element  formulation 
and  then  specialized  to  a  cylindrical  shell  geometry. 
Linear  results  are  compared  to  exact  elasticity  solutions  of 
a  laminated  cylindrical  shell  in  cylindrical  bending.  The 
approach  predicts  displacement  and  stress  for  the 
unidirectional  cases  very  we 1 1  even  for  the  very  thick 
laminates . ^.Transverse  stress  for  the  cross  ply  cases  are 
not  as  accurate  Since  equilibrium  at  the  ply  interfaces  is 
not  enforced,  a  typical  drawback  of  2-D  theories. 
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effects.  As  the  shell  becomes  thicker  relative  to  its  radii 
of  curvature  and  in-plane  dimensions,  or  for  laminated 
composite  shells,  these  transverse  effects  become  more 
pronounced,  especially  the  transverse  shearing.  The  first 
theories  that  included  the  transverse  shear  deformation 
relax  the  assumption  on  the  deformed  normals  of  the  shell 
middle  surface.  The  shell  kinematics  introduced  by  Bassett, 
as  discussed  in  [11],  express  the  shell  displacements  as  an 
infinite  power  series  in  the  thickness  coordinate. 
Hildebrand,  Heissner,  and  Thompson  C 1 1 3  introduced  truncated 
Bassett  kinematics  to  analyze  the  importance  of  the 
transverse  stresses  and  strains  in  thin  elastic  orthotropic 
shells.  Maghdi  [12]  applied  similar  truncated  series 
representations  for  general  thin  isotropic  elastic  shells. 
Hildebrand,  et  al ,  found  that  the  effects  of  the  second 
order  displacement  terms  and  terms  in  the  transverse 
displacement  that  give  a  nonzero  transverse  normal  strain 
were  negligible.  Heissner  used  similar  kinematics  to 
analyze  plates  [13]  and  then  sandwich  shells  [14].  Mindlin 
similarly  included  rotatory  inertia  terms  in  the  dynamic 
analysis  of  plates  [15].  The  preceding  studies  have 
resulted  in  the  so  called  Beissner-ifindl in  CRM)  theory.  The 
approach  does  not  satisfy  the  transverse  shear  boundary 
conditions  on  the  top  and  bottom  surfaces  of  the  shell  or 
plate  since  a  constant  shear  angle  through  the  thickness  is 
assuned ,  i.e.,  plane  sections  remain  plane.  Because  of 
this,  the  theories  based  on  these  kinematics  usually  require 
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shear  correction  lactors  for  equilibrium  conajosnoiona. 

iatij  .arris,  and  Stavsky  Lid  generalized 
Beissner-kiindlln  piaca  theory  to  laminated  composite  plates. 
Whitney  and  Paganc  117J  were  fi ral  to  app ^ .  i  to 
composite  plate  analysis.  Beddy  1 1 3  J  n as  appiieo  hi,  theory 
to  linear  anisotropic  shell  structures  ox  constant  principal 
and  twist  curvatures.  door  il&j  applies  Rt»  tusory  to 
examine  the  Stability  of  laminated  plates.  Thick  composite 
plate  closed  form  solutions  are  reported  by  Beddy  and  Chao 
[20].  Linear  static,  stability,  and  vibi ation  analyses  of 
laminated  plates  and  shells  were  pencr-med  by  Boor  and 
Mathers  [21].  Two  finite  element  i ormulat i ons ,  mixed  and 

displacement,  are  compared.  Ahmad ,  Irons,  and  Zienkiericz 

i 

[221  introduced  a  finite  element  approach  with  independent 
transverse  displacement  and  rotational  degrees  of  freedom 
such  that  a  Beissner -Hi  nol  iv\  shear  deformable  shell  element 
results.  These  elements  perform  admirably  for  thick 
situations,  tut  develop  shear  locking  numerical  difficulties 
fc-  thin  cases  and  consequently,  reduced  elemental 
Integrations  [23]  or  tlscrev«  Kirchhofl  constraints  [24] 
become  necessary . 

Levinson  [25],  Wurth}  126],  and  Beddy  [27]  have 
developed  plate  theories  that  include  cubic  terms  in  the 
in-plane  displacement  kinematics.  Satisfying  zero 
transverse  shear  stitt-  on  uor  top  anu  bottom  surfaces  of 
the  plate  results  in  a  parabolic  shear  strain  distribution 
through  the  thickness.  t.h  agreeing  more  closely  with 
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elasticity  results. 


The  number  of  variables  in 


tbe 


kinematics  is  equal  to  that  in  tbe  RM  theory,  but  shear 
correction  factors  are  not  required.  Reddy  and  Phan  [28,29] 
analyzed  linear  laminated  anisotropic  plates  using  the 
parabolic  shear  theory.  The  parabolic  shear  theory  is 
extended  to  shell  geometries  in  a  study  of  the  linear 
behavior  of  laminated  shells  by  Reddy  and  Liu  [301. 
Bhimaraddi  has  applied  the  parabolic  shear  strain 
distribution  to  analyze  the  linear  vibrational  behavior  of 
isotropic  cylindrical  shells  [31]  and  rings  and  arches  in 
[321.  Additionally,  Soldatos  applied  the  parabolic  shear 
theory  to  examine  the  stability  and  vibrations  of  laminated 
circular  [33]  and  noncircular  [34]  cylindrical  shells. 

Geometrically  nonlinear  behavior  is  captured  by 
allowing  large  displacements  and  rotations  in  the  shell 
strain  displacement  relations.  An  approach  that  neglects 
the  nonlinear  terms  that  represent  in-plane  rotations  of  the 
shell  is  an  example  of  the  intermediate  nonlinear  theories 
often  used  in  stability  analysis  [35,36].  Small  strain, 
moderate  rotation  approaches  for  anisotropic  shells  are 
given  by  Librescu  and  Schmidt  [37]  and  Schmidt  and  Reddy 
[70].  Successive  approximations  made  to  exact  shell  strain 
displacement  relations  where  large  strains,  displacements, 
and  rotations  are  all  initially  allowed,  are  presented  for 
isotropic  shells  by  Sanders  [38]  and  Pietraszkiewicz  [71] 
and  anisotropic  shells  by  Librescu  [39]. 

An  '*arly  work  by  °chtnit  and  Monf orton  [40J  formulates 

% 


an  anisotropic  cylindrical  shell 


element  based 


upon  the 


classical  theory  that  allows  intermediate  geometric 
nonlinearities.  Stolarski  et  al  [41]  present  a  simple 
triangular  shell  element  formulation  that  includes 
interswdiate  nonlinearity.  Beddy  and  Chandrashekhara  [42] 
solve  both  cylindrical  and  spherical  laminated  shell  cases 
assuming  RM  transverse  shear  and  an  intermediate 
nonlinearity.  Beddy  [43]  analyzed  the  nonlinear  transient 
behavior  of  composite  plates  including  RM  transverse  shear. 
Futcha  and  Reddy  [44]  used  the  parabolic  shear  theory  in 
formulating  a  mixed  element  for  nonlinear  anisotropic  plate 
analysis.  Moor  and  Peters  [45]  give  the  nonlinear  response 
of  anisotropic  cylindrical  panels  via  a  Hu-Washizu  mixed 
shallow  shell  finite  element  approach  that  includes 
transverse  shear  deformation  together  with  a  Rayleigh-Ritz 
technique.  Meroueh  1 4.0]  e.r.o  Scrana  [47,48]  develop 
geometr leal ly  nonlinear  shell  finite  element  approaches  that 
are  shear  deformable.  Stein  [49]  uses  truncated  series 
expansions  of  exact  nonlinear  strains  in  a  shell  approach 
that  also  includes  transverse  shear  deformation. 
Qeosietrical ly  nonlinear  quasi  •  three -di mens i  onal  approaches 
for  laminated  composite  plates  and  shells  have  been 
developed  by  Palazotto  et  al  [50.5;]. 

Very  recent  laminated  plate  and  shell  review  articles 
are,  respectively,  due  to  Bo or  and  Burton  [68]  and  kapania 
[69] . 


The  present  report  presents  an  approach  for  laminated 
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composite  shell  analysis  that  combines  a  parabolic 
transverse  shear  distribution  with  a  simplified  large 
displacement  and  rotational  geometric  nonlinearity. 
Solutions  to  the  theory  are  found  via  a  displacement  based 
finite  element  approach.  Extensive  results  are  given  for 
the  linear  analysis  of  a  laminated  cylindrical  shell  in 
cylindrical  bending.  This  simple  case  is  chosen  to 
illustrate  the  effectiveness  of  the  2-D  approach  in 
describing  3-D  response.  Once  its  limitations  are  known, 
the  approach  can  then  be  applied  to  the  majority  of 
situations  where  elasticity  solutions  are  not  possible. 


1 1 .  Theoretical  Development 


General 

Although  shell  theory  can  be  based  entirely  on  surface 
definitions,  transverse  shear  strains  and  stresses  are  then 
not  easily  included.  One  way  to  define  strain  displacement 
relations  that  can  easily  incorporate  three  dimensionality 
(3-D)  is  to  specialize  the  general  3-l>  strain  displacement 
relations  expressed  in  arbitrary  orthogonal  curvilinear 
coordinates.  To  this  end,  consider  the  infinitesimal  line 
segment,  MH ,  of  length  ds  embedded  in  a  differential  volume 
element  in  Figure  1.  This  differential  volume  is  linearly 
transformed,  i.e.,  deformed,  to  a  new  configuration  where 

ft 

the  line  segment  is  now  of  length  ds  .  As  a  result  of 

*  * 

deformation,  the  line  segment  MN  of  Figure  1  moves  to  MM 
represented  by  the  displacement  vector,  u.  By  subtracting 
the  original  and  deformed  squared  lengths  of  the  line 
segment,  the  Green's  strain  tensor,  j  (i,j=1.2.3),  is 
defined  as  shown  in  Eqn  (1) . 


(ds  ) ^  -  (ds)  ^  =  1y  .  .  dy . dy  . 

i  J  l  J 


(repealed  indices  sum) 


(1) 


where  the  y are  the  orthogonal  curvilinear  coordinates  of 
the  undeformed  system.  The  physical  strains,  e  ,  are  then 
found  from. 


e 


i  j 


yiJ 


h  .  h  . 
i  J 


(no  sum) 


(2) 


In  Eqn  ),  the  b^  are  scale  factors  defined  from  the 
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metric  tensor  of  the  undeformed  coordinate  system  and  the 
are  shown  in  the  Appendix,  Eqn  (Al) ,  where  the  ir  are 
the  coordinates  of  the  displacement  vector,  u.  For  the  case 
of  small  strains  (£  <  .04),  the  c (no  sum)  are  the 

elongations  of  the  fibers  of  the  differential  volume  element 
and  the  (i^j)  are  the  shears,  i.e.,  the  difference  from 

ninety  degrees  originally  perpendicular  fibers  are  oriented 
after  def ormation .  The  elongations  and  shears  are 
identically  zero  for  rigid  body  movements  of  the  element 
[52]  . 


FIGURE  1.  Segment  MN  Deforms  to  11*11*  Through  Displacement 
Vector  u. 

Now  consider  a  shell  geometry  that  can  be  described  by 
orthogonal  curvilinear  middle  surface  coordinates,  £ ^  and 


*0,  surface  normal  C.  and  radii  of  curvature,  Rj  and  R^  as 
shown  in  Fi fore  2.  For  this  geometry,  the  ncale  factor 
terms  of  Eqn  (2)  are  given  by, 


\  =  V*  -  C'V’  h2  =  V* 

where  a  (y=l,2)  are  related  to  the 
Y 

metric  l *  1  , 


C/Ha) ,  h3  =  1  (3) 

cements  of  the  surface 


FIQUHK  2.  Arbitrary  Shell  Coordinate  System. 

For  a  body  of  volume  V  in  equilibrium  with  prescribed 
forces  on  a  part  of  its  surface,  S, ,  and  prescribed 
geometric  boundary  conditions  on  the  remaining  surface,  S^. 
that  has  undergone  an  infinitesimal  virtual  displacement, 
6u ,  we  have , 
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P,  6u,  )  dV 

k  k 


F.  6u,  dS  =  0 

k  k 


(4) 


|  ‘"u^U 


'I 


where,  a.  .  are  elements  of  the  2nd  Piola-Kirchhof f  stress 
i  J 

tensor,  c .  .  are  defined  in  Eqn  (2) ,  P  are  components  of 

1  J  ^ 

prescribed  body  forces,  and  are  components  of  prescribed 
surface  tractions  [53].  Assuming  the  existence  of  a  strain 
energy  density  function.  W* ,  small  strain  deformation,  and  a 
linear  constitutive  law,  the  principle  of  stationary 
potential  energy  can  be  stated  as, 

6n  =  0  (5) 

p 


and  if  the  forces  vary  neither  in  magnitude  nor  direction 
during  virtual  displacements, 


n  =  J  W*  ( u .  )  dV  -  f  F.  u.dS  (6) 

P  JV  1  JSj  1  1 

ft 

where  W  (u.)  =  1/2  a.  ,  a..,,  are  elastic  constants, 

l  ijkl  ij  kl  ljkl 

In  Eqns  (5,6),  the  strain  energy  density  is  written  in 
terms  of  the  displacement  components  by  using  Eqns  (2,A1). 
The  integrals  are  taken  over  the  original ,  undeformed  volume 
and  surfaces  consistent  with  a  Lagrangian  viewpoint. 
Additionally,  the  array  of  elasticity  constants,  are 

defined  with  respect  to  the  original  configuration.  T53J 

All  of  the  necessary  tools  have  been  introduced  for  a 
Lagrangian  displacement  based  energy  approach  to  solve 
elastic  geometrically  nonlinear  shell  problems.  To  develop 
shell  strain  displacement  relations,  one  should  first 
identify,  respectively  the  general  3D  curvilinear 


coordinates  y^(  y? ,  and  y^  of  Eqns  (1.2.A1)  with  the  Shell 
coordinates.  and  £  of  Figure  2.  Then,  by  defining 
kinematics  that  rexate  the  middle  surface  displacements  to 
the  continuum  displacements  of  the  shell.  strain 
displacement  relations  are  written,  and  finally,  the  total 
potential  energy  expressed  in  terms  of  displacement  using 
Eqns  (5.6)  is  extremized  with  respect  to  displacement 
resulting  in  the  equilibrium  equations.  The  definition  of 
the  kinematics  for  these  types  of  formulations  is  obviously 
very  important.  The  assumptions  that  will  eventually  lead 
to  general  shell  kinematics  and  a  potential  energy 
expression  are  next  discussed. 

Basic  Assumptions 

A  geometrically  nonlinear  theory  governing  the  shell  is 
based  on  the  following  assumptions: 

1.  The  shell  is  thin  and  therefore  assume  that  it  is  in  an 
approximate  state  of  plane  stress,  i.e.,  the  transverse 
normal  stress  is  approximately  equal  to  zero.  This 
assumption  effectively  reduces  the  generally  three 
dimensional  behavior  of  the  shell  so  that  it  can  be 
described  by  the  behavior  of  only  a  datum  surface;  this 
results  In  a  2-D  approach. 

2.  The  transverse  shear  stress  distribution  is  parabolic 
through  the  shell  thickness  and  vanishes  on  the  top  and 
bottom  surfaces  of  the  shell. 

Because  the  shell  is  thin..  in -plane  stress  tends  to 
dominate  the  shell  response  under  loading,  i.e.,  the 
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transverse  stresses  are  of  lesser  importance.  However,  it 
is  well  known  that  laminated  thin  flat  plate  response  can  be 
significantly  influenced  by  transverse  shear  deformation. 
Similar  results  have  been  found  for  laminated  shells 


[30,541.  Therefore,  transverse  shear  stress  is  included  in 
the  present  model,  via  a  2-D  approach  whereas,  transverse 
normal  stress  is  not.  This  approach  is  consistent  with 
remarks  made  by  Koiter  [10,55]  and  John  [561  who  performed 
order  of  magnitude  studies  on  the  stresses  in  general  shell 
geometries.  Additionally,  their  results  are  seen  to  hold 
generally  true  in  the  elasticity  solutions  of  both  laminated 
flat  plates  by  Pagano  [57,58]  and  laminated  shells  in 
cylindrical  bending  by  Ren  [59]. 

3.  The  shell  is  restricted  to  small  strains  and  consists  of 
linear  elastic  laminated  orthotropic  material. 

4.  Exact  Green's  strain  displacement  relations  are  assumed 
for  the  in-plane  strains  ^22'  and  *12  of  E<lns  (2.3)  and 
linear  strain  displacement  relations  are  assumed  for  the 
transverse  strains  £:  and  e  . 

■o  vJ  lu 

5.  Since  the  shell  is  assumed  thin,  the  normals  to  the 
middle  surface  are  approximately  inext.ensibl  e ,  hence  the 
transverse  displacement  is  constant  through  the  thickness. 
Consistent  with  the  preceding  assumption,  retention  of  only 
linear  strain  displacement  terms  then  gives  -  0. 

Librescu  [39],  as  mentioned.  develops  general  shell 
equations  based  upon  varying  approximations  in  the 
magnitudes  of  both  strain  and  the  rotations  differential 
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elements  undergo  during  deformation.  A  consistent  small 
strain,  moderate  rotation  theory  has  nonlinear  (not  exact!) 
strain  displacement  relations  for  the  in-plane  strains  yet 
linear  relations  for  the  transverse  shear  strains.  The 
present  approach  can  therefore  be  viewed  as  a  small  strain, 
simplified  large  rotation  theory  due  to  the  exact  in-plane 
strain  assumption.  The  accuracy  in  rotation  is  limited  by 
the  linear  (in  displacement)  assumption  on  the  transverse 
shear  strains. 

Constitutive  gelations 

Materials  consisting  of  unidirectional  fibers  embedded 
in  a  matrix  are  of  interest  in  this  study.  For  fibers 
alligned  with  the  1  material  axis,  we  can  assume  that  the 
material  is  transversely  isotropic  (engineering  constants 
E2=E3  and  yj2=v>i3^  with  respect  to  planes  parallel  to  the 
2-3  plane  [60],  In  this  case,  the  constitutive  relation 
used  in  Eqn  (6)  becomes , 


C11 

°12 

°13 

0 

0 

0 

°2 

C12 

C22 

C23 

0 

0 

0 

*2 

*3 

C13 

C23 

C33 

0 

0 

0 

*3 

r  S 

0 

0 

0 

°44 

0 

0 

C4 

0 

0 

0 

0 

C55 

0 

*5 

>. 

0 

0 

0 

0 

0 

cee 

(7) 


where  contracted  notation  is  introduced,  i.e.,  <?  ^  -  (no 


sum) ,  «. 


°23 ’  °5  =  ‘’’IS  °6  =  *12'  and  *i  =  *ii 


(no  sum)  , 


c .  =  2c..,  c. 
4  23  5 


2£13-  *=6  =  2<r  1 2 


Using  assumption  1  (o^  =  0)  in  Eqn  (7),  solve  for 
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1 


(R) 


Substituting  back  into  Eqn  (7) ,  £  is  eliminated,  and  the 
constitutive  relations  in  material  axes  become. 


Ci3Cj3 


where  Q . .  =  C . .  - 

i  J  i  J 


.  and  for  transverse  isotropy,  it 


is  easily  shown  that,  Qjj-Ej/A.  Q12=U21E2/A>  Q22=E2/A’ 

Q66=:Q12’  Q44=G23’  Q55=Q13'  A  1  “  V12V21' 

Consider  a  shell  that  is  constructed  of  layers  of  the 


transversely  isotropic  material  described  by  Eqn  (9). 
Generally,  the  fibers  of  the  ktb  individual  layer  or  ply  are 
oriented  at  an  angle  with  respect  to  the  middle  surface 
shell  coordinates  and  Therefore,  tha  constitutive 
relations  of  Eqn  (9)  for  that  ply  must  be  transformed  into 


shell  coordinates  resulting  in  Eqn  (101. 


N 

k  r 

h 

'  35 

Q11  ®12  ®] 
°22  ®5 


'66 


O'-  r~  3‘  ti;l 


k 

r*.i 

-4 

k1 

CIO) 


•whfi7‘&  ,  the  Q,  .  (i,j -1,2.0)  and  Q  (a,n=i,5!  are  elements  of 
i  j  mn 

symmetric  arrays  of  transformed  stiffnesses  for  the  kth  ply 
and  a .  ,  c  ,  a  ,  and  £  are  measured  with  respect  to  shell 

X  1  fll  IB 

coordinates  Z  and  C> 

a 


K  j  C3 

Consider  the  following  kinei^&tics  for  the  arbitrary 
shell  described  by  orthogonal  curvilinear  coordinates, 
referring  back  to  Figure  2  for  conventions. 


UjCCj.fj.O  =  U(1 

-  C/B  > 

+  r«/1 

r3 

+  C  J'l 

♦  c4* 

U2(W°  =  V(1 

'  C/B2) 

+  <>2 

* 

,  r  3 
*  *  r7. 

♦  cV 

U3(W  =  W  (11> 

where,  u,  v,  w ,  V  ,  <P  ,  Y  ,  and  6  ,  are  functions  of  and 

a  a  'a  a  1 

v  are  rotations  of  the  normals  and;  <£  ,  r  .  0  ,  are 

i  a  a  a  a 

functions  to  be  determined  such  that  and  vanish  on  the 
top  and  bottom  surfaces  of  the  shell. 

The  assumed  displacements  of  Eqn  (11)  are  slightly  more 
involved  than  those  introduced  for  flat  plates  by  [271.  For 
a  plate,  only  third  order  terms  in  the  thickness  parameter. 
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were  necessary  to  give  the  desired  parabolic  transverse 
stress  distribution.  Shell  structures,  due  to  their  curved 


surfaces,  however,  have  coupling  between  displacements  that 
plates  do  not  and  the  fourth  order  terms  become  necessary. 
The  kinematics  will  ultimately  closely  resemble  the  plate 
kinematics  that  have  been  used  in  the  literature. 
Additionally,  Reddy  [30]  arrives  at  identical  general  shell 
kinematics  as  derived  here  but  uses  Sanders  shell  relations. 

Keeping  only  linear  displacement  terms  from  Eqns 
( 2 , A 1 ) ,  the  transverse  shear  strains  become. 


(u  o  +  h  U  •*  “  uoho 

U  ,  A  4  A  ,  «J  /  <6  ,  U 


-t—  (u,  ,  ♦  h.u.  _  -  u.h.  _) 
u  j  3,1  I  1,3  Ll.,3 


(12) 


■her#  (  ) ,  (ot=l,2)  refers  to  differentiation  with  respect 

to  ?  and  (  ),,  refers  to  differentiation  with  respect  to  £ . 

Of  o 

Using  Eqn  (11)  in  the  first  of  Eqn  (12)  and  assuming 
zero  transverse  shear  stress,  and  therefore,  strain,  on  the 
Shell  boundaries  (see  Eqns  (9,30))  the  following  is  derived. 


*2  =  °- 


e  =  — 

2  2R„  ’ 


1  - 


8R 


2-> 


♦  -LI  ) 

3h2  2  “a 


where  h  is  the  shell  thickness,  see  Figure  2. 

Assume,  for  the  moment,  that  we  have  a  fairly 
shell,  i  .  ,  let  R  -  5h  In  thin  case  the  underlined 


(13) 


thick 

term 


in  Eqn  (13)  is  only  .005  and  therefore  can  be  neglected 
compared  t«,  1  as  indicated.  Furthermore,  we  can  neglect  the 
fourth  order  term  (#9)  of  Eqns  (11,13)  since  it  is  only  1/20 
of  the  third  order  term.  A  similar  exercise  is  applied  to 
£5  giving  for  the  general  shell  kinematics. 


vw° 

u2(?1.?2.C) 


C  3  W'l 

-57  1  *  Oi  *  <  ‘""i  *  -57  > 


V*  -  T  1  *  <>2  *  ^^2  *  TT  > 

4  * 


VW  =  " 


(14) 


where ,  k  =  — . 

3h 

The  rotation  of  the  normal.  V'  .  is  a  result  of  bending 

c* 

deformation  and  therefore,  the  slope  of  the  elastic  curve, 
w.  ^ ,  is  always  a  combination  of  y  and  shear,  as  shown 

in  Figure  3.  For  a  flat  plate,  C/B  =0  and  a  =  1  and  Eqn 

a  y 

(14)  reduces  to  that  used  in  the  literature. 
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She! 1  Strain  Displacement  Relations 

The  kj.-.nmatics  of  Eqn  (14)  next  define  the  general 
shell  strain  displacement  relations  using  Eqn  (Al)  as  given 
in  Eqn  ( 15)  . 

C .  =  £°  *  Cp  *.  ,  i  =  1 ,2.6  (15) 

1  1  i  p 

p  dun  1  to  7 

*.  =  £m  +  ^  *m2’  “=4*5 

where,  e°  x, ,  through  £°  ,  x  „  are  shown  in  the 

Appendix,  Eqns  (A2-A6) ,  specialized  for  cylindrical 
coordinates.  Mote  that  the  convention  on  i  (and  j)  has  been 
changed  so  that  instead  of  taking  on  values  1,2,  or  3,  they 
now  represent  the  values  1.2,  and  6.  This  will  always  be 
the  case  from  this  point  on.  The  first  subscript  of  x 


refers  to  the 

strain  component  I,*! 

.  or ,  8  and 

the 

second 

subscript  refers  to  the 

power  of  r 

it  multiplies. 

The  e ‘ 

and  the  x^ 

through 

terms 
i  i 

consist  of 

displacement 

expressions 

that  are 

dependent 

only  on 

the 

surface 

parameters  ?  ^ 

and  ?2. 

In  developing  the  2-D  theory,  it  is  convenient  to 
define  the  strains  such  that  the  strain  terms  multiplying 
the  C ' s  in  the  expressions  of  Eqn  (15)  are  functions  only  of 
the  in-plane  parameters,  ( ^  and  (j.  Therefore,  in  arriving 
at  Eqn  (15),  binomial  series  expansions  such  as  those  shown 
in  Eqn  (16)  were  used  for  this  purpose.  Terms  similar  to 
the  left  band  sides  of  Eqn  (16)  result  dtie  to  the  shell 
scale  factors,  i.e.  those  terms  in  Eqns  (A2-A4)  that  are 
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1 - ■M*  !■» 


t  7T - - 


functions  of  h..  For  the  thick  shell  case  where  B  =5h  and 

l  ct 

C=h/2,  the  first  two  terms  on  the  right  hand  sides  for  each 
expression  in  Eqn  (16)  represent  99X  and  97X ,  respectively, 
of  the  total  on  the  left  hand  sides.  Consequently,  series 
expansions  as  in  Eqn  (16)  and  similar  expansions  are 
truncated  after  the  first  order  C  terms. 


1 


1 


H  (1-C/R  ) 

a  a 


(16) 


r  i 

A, 

1 

♦  rem 

1  R, 

B„ 

*•  1 

2 

J 

For  simplicity,  in  arriving  at  the  transverse  shear 
strains  of  Eqn  (15) ,  the  cubic  term  found  as  a  result  of 
Eqns  (12,14)  is  neglected  as  small  using  the  R^  =  5h 
argument  as  before.  This  results  in  a  symmetric  transverse 
shear  strain  and  stress  distribution  about  the  midplane  of 
the  shell.  In  reality  however,  as  the  shell  becomes 

thicker,  the  transverse  shear  stress  tends  to  become  less 
symmetric  about  the  midplane. 

Shell  Potential  Energy 
From  Eqn  (6) ,  let 

n  =  o  +  v  (17) 

p 


where,  U  is  the  internal  3train  energy,  and  V  is  the  work 


70 


done  by  the  external  forces. 

Bqns  '10,17,1b)  then  give  for  the  internal  strain 
energy  of  the  arbitrary  shell. 


~*\Ci 

ii 

eH 

D 

u 

[9nu! 

*  fp«.P>2  *  4 

■f 

a5ia'*I 

*  V  <*2  *  *  *«.'«;  *  fP)V2 

*«!«<*; 

*  '"‘ip'K  *  <'■„> 

♦ 

2^26U2 

*  4  '’‘’'Or’]  dC<10  • 

(18) 


°2  = 


*JJ  f 


044(iC4  +  C  *42J  *  Q55U5  +  C  *52J 


O  h 


*  2Q45(<  ♦  -  <%2>]  dCdO 

where,  D  =  U  ♦  0,  p,r  -  1,2, 3, 4, 5, 6, 7,  Q  generally  vary 
1  2  I*/ 

as  a  function  of  C  since  a  laminate  is  constructed  of  plies 
with  various  fiber  orientations,  and,  0  represents  tbe  shell 
middle  surface  - 

Manipulations  of  Frjn  (18)  are  facilitated  if  the 
following  quantities  are  defined.  From  Eqn  (15),  let 


*  =  *°  +  7C  Z 


(10) 


where , 
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Using  Eqn  (10)  in  the  first  of  Eqn  (18),  we  get. 


-  -5-  J  |  10  *  >T  * 


dCdO 


=  -i—  [  [  -[  c‘  c°  Q  +  2  c"  Q  x  CP 

2  J  J  \  j  i  wij  j  wij  ip 

n  h 

*  “ip  *ir  «U  fP*r  «  }  «" 

Q11  fl2  ®16 

where,  <2  =  Q22  Q2g  ,  i,j  =  1,2,6  and  p,r  sub  1  to  7 . 

Q66J  superscript  t  =  transpose 

The  second  integral  of  Eqn  (20)  bints  at  the  complexity  of 
the  in-plane  strain  energy  of  the  shell.  The  energy  is 
comprised  of  many  strain  squared  terms  due  to  summations  on 
up  to  four  different  indices.  Finally,  Eqn  (20)  is 
concisely  written  in  terms  of  an  area  integral  representing 
the  shell  midsurface  as  shown  in  Eqn  (21)  where  the  C 
dependence  has  been  integrated  by  defining  a  series  of 
elasticity  arrays  shown  in  Eqn  (22).  The  of  Eqn  (21)  are 
written  in  terms  of  strain  components  and  elements  of  the 
elasticity  arrays  in  the  Appendix,  Eqn  (A7) . 


?2 
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(21) 


0 


1 


+  ^2  +  U3  }  d° 


CAlj’  Bij’  Dij*  Elj'  Fij’  °ij-  Hij'  Tij*  Jij*  Kij' 


I 


L. . ,  P. . ,  H. . ,  S. . ,  T. . ]  = 

U  ij’  lj’  i j  u 


QAj  Cl.C.C2^3^4^5^6^7^8^9^10^11^12^13^14]  dC 


(22) 


•here  i.j  =  1,2,  and  6. 

A  similar'  manipulation  is  performed  on  the  second  of  Eqn 
(16),  the  shear  terms,  giving. 


°J  =  T-  J  1 


(  e°£°  A  ♦  2c° a  _  D  +  a  _  F  )  dft 
mnmn  n  m2  mn  n2m2mn 


and,  with  m,n  =  4,5,  (23) 

t A  . D  . F  ]  =  T  0  [  l.C2.C4)  dC  (24) 

mn  mn  mn  J  on 

h 

Primarily  due  to  the  assumed  nonlinearity,  many  new 
elasticity  arrays  are  shown  above.  The  energy  terms  that 

contain  the  higher  order  elasticity  constants  are  also 

2 

multiplied  by  powers  of  k=-4/3h  (see  Eqn  (14)).  Because  of 
the  squared  shell  thickness  term  in  the  denominator  of  k, 
the  energy  terms  are  not  as  many  orders  of  magnitude  apart 
as  Eqn  (22)  may  at  first  indicate.  For  example,  the  , 

F  ,  and  energy  terms  are  all  the  same  order  of 

magnitude  in  h  or  (  despite  4  orders  of  magnitude  difference 
in  the  definitions  of  the  elasticity  arrays  given  in  Eqn 
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T 


7^ 


(22)  . 

In  a  general  geometrically  nonlinear  analysis,  the 
energy  due  to  the  externally  applied  loads  contains 
nonlinear  displacement  terms.  These  arise  when,  because  of 
large  rotations,  the  applied  loading  develops  higher  order 
coupling  components  in  coordinate  directions  where  no 
loading  originally  existed.  As  discussed  by  Brush  and 
Almroth  [35],  these  higher  order  terms  give  potential  energy 
contributions  that  are  negligibly  small  for  the  cases  of 
intermediate  nonlinearity,  i.e.,  moderate  rotations. 
However,  the  large  rotational  case  generally  must  account 
for  them.  If  the  loading  is  actually  a  prescription  of 
displacement,  the  higher  order  nonlinear  loading  terms  need 
not  be  included  since  applied  loading  is  zero.  This 
alternative  simpler  approach  is  taken  here. 

In  summary,  Eqn  (21,23)  with  Eqns  (A7 , 22 , 24 , 15 , A2-A6) , 

give  the  total  potential  energy  of  the  arbitrary  shell, 

where  a  specific  Shell  geometry  is  defined  by  Specifying  the 

scale  factor  terms  of  Eqn  (3) .  This  expression  is  then 

extremized  with  respect  to  the  displacements  giving  the 

equilibrium  equations.  The  result  of  the  extremization  is 

five  coupled  nonlinear  partial  differential  equations.  The 

present  approach  gives  nonlinear  extensional  terms,  i.e., 

terms  in  A. .,  as  well  as  nonlinear  terms  in  D. .  through  T. . 

ij  lj  °  ij 

due  to  the  nonlinearity  in  the  curvature  terms  (x^)  of  Eqn 
(15).  The  moderate  rotational  theories  retain  in  the  strain 
displacement  relations  nonlinear  displacement  terms  only  in 
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**  "T 


Tjr 


T 


the  transverse  displacement,  w ,  as  a  result  of  assuming  that 

the  rotations  are  small  compared  to  unity  [361  This 

results  in  nonlinearity  in  the  equilibrium  equations  in  the 

extensional  A.,  terms  only. 

1  j  ’ 
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III. 


Finite  Element  Solution 


The  variation  of  the  potential  energy,  n  ,  of  the 
discretized  domain  gives  coupled  nonlinear  algebraic 
equations  where  the  unknowns  are  the  values  of  displacement 
at  the  nodes.  In  a  generally  nonlinear  analysis,  these 
equations  are  not  usually  solved  directi v  but  instead,  are 
linearized  and  solved  by  incremental/iterative  methods.  The 
linearized  equations  are  found  by  effectively  carrying  out 


an  additional 

differentiation  of 

n 

p 

with 

respect  to  the 

displacements . 

x'he 

man i pul at ions 

that 

follow  give  a 

convenient  form  for 

0  such 

p 

that 

the 

variation  and 

1 inear ization 

are 

straight 

forward 

to  accomplish. 

Afterward,  a  36  degree  of  freedom  (dof)  curved  element  is 
discussed . 

Element  Independent  Formulation 

We  seek  a  form  for  the  potential  energy,  n  ,  shown  in 

P 

Eqn  (2S)  ,  such  that  the  first  variation,  <511^,  is  given  by 
Eqn  (26)  and  the  linearized  incremental/iterative  equations 
are  given  in  the  last  of  Eqn  (28) .  Eqn  (28)  is  found  from 
Eqn  (26)  by  expansion  of  the  equilibrium  equations,  F(q), 
into  a  Taylor  aeries  about  a  small  Aq  and  truncating  as 
shown  in  Eqn  (27) . 


n 

p 


T 

q  b 


(25) 


6qT  F  (q)  =0  (26) 
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where,  q  is  a  column  array  of  nodal  displacements,  B  is  a 
column  array  of  nodal  loads,  K  is  an  array  of  constant 
Stiffness  coefficients,  lf^  is  an  array  of  stiffness 
coefficients  that  are  linear  in  displacement,  is  an  array 
of  stiffness  coefficients  that  are  quadratic  in 
displacement . 

From  Eqn  (26),  for  arbitrary  and  independent  6q,  F(q)  = 
0 ,  then , 

d  F 

F(q  +  Aq)  =  F(q)  +  -  Aq  ♦  ...  =0  (27) 

d  q 

9  F 

-  Aq  =  [  K  ♦  +  JJz  ]  Aq  =  KrAq  =  -F(q)  (28) 

a  q 

•here  Kt  is  the  tangent  stiffness  matrix. 

Eqn  (28)  gives  solutions  in  an  iterative  manner  via  a 
Bewton-Baphson  technique.  The  current  values  of  the 
elements  of  q  are  substituted  into  ■  and  V  resulting  in  an 

i  2 

array  of  constants  for  Kt<  The  set  of  linear  equations  are 

then  Solved  for  Aq,  which  is  then  added  to  q  giving  the 

updated  nodal  displacements.  Assuming  the  solution  is  not 

yet  converged,  the  right  hand  side  (RHS)  of  Eqn  (28)  is 

nonzero.  Iteration  continues  until  F(q)  becomes  arbitrarily 

small  signifying  equilibrium  has  been  satisfied. 

The  preceding  general  development  assumed  a  certain 

formalism  in  the  definitions  of  the  arrays,  K,  V  ,  and  N  in 

1  2 

np  of  Eqn  (25)  such  that  these  terms  repeat  themselves  in 
the  first  variation  of  (Eqn  (26))  and  again  in  the  Taylor 
ser ies  expansion  that  gave  the  incremental  equations  (Eqn 
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(28)).  Ensuring  this  formalism  mill  greatly  decrease 


the 


required  number  of  algebraic  manipulations  since  once  K,  H  , 
and  Vz  are  formed  in  the  expression  of  Eqn  (25)  ,  the 
variation  and  linearization  processes  become  straight 
forward.  Additionally,  fewer  programming  lines  are  needed 
since  the  same  arrays  are  present  in  both  the  equilibrium 
and  incremental  equations  and  hence  need  be  calculated  only 
once  per  iteration  in  solving  Eqn  (28)  . 

Rajasekaran  and  Murray  [61]  outline  a  procedure  that 
will  give  the  desired  repetition  and  also  symmetry  of  the 
three  components ,  K,  If  ,  and  N  ,  that  comprise  K  .  Some 

12  T 

generalizations  are  used  in  the  present  case.  We  start  by 
dividing  each  strain  component  of  Eqn  (15)  into  linear  and 
nonlinear  parts  as  shown  in  Eqn  (29) . 


c  =  l;  d  ♦ 

i  o  i 


*.  =  L,  d  + 

ip  pi 


1  ,T  _ 

~n —  d  _  H .  d 

2  o  i 


4-  dT  H.  d 

2  pi 


(29) 


"here.  L.  (J  =  0,1, 2,... ,7)  are  column  arrays,  H. 

w  1  J  1 


ape 


symmetric  matrices,  d  is  the  displacement  gradient  vector, 
i  =  1,2,6  and  p=l  to  7,  as  before,  for  p=5,6,7  contain 
all  zeroes. 

The  displacement  gradient  vector  contains  the  unique 
displacement  terms  that  define  the  strain  components  in  Eqn 
(IS).  Since  there  are  18  unique  displacement  terms  in  Eqn 
(15),  the  vector,  d,  is  an  18x1  array  as  shown  in  Eqn  (30). 
The  way  in  which  the  elements  of  d  are  approximated  defines 
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the  specific  element  type,  i.e.,  number  of  dof ,  shape 

functions,  etc.  This  elemental  definition  is  done 

independent  of  the  repeating  formalism  of  Eqns  (25,26,28). 

The  L,  and  H.  therefore,  are,  18x1  and  18x18  arrays 
j  i  J  i 

respectively,  and  contain  constants. 


iT=  ju  u ,  j  u ,  2  v  v  ,  l  v , , 


' 1  *2  *  11  *22  ’  12 


*1  *1.1  *1.2  *2  *2.1  *2.2  } 


(30) 


Next,  Eqn  (30)  is 

substituted 

into 

the 

internal 

strain 

energy  of  Eqn 

(21)  giving 

for 

the 

in-plane 

energy 

expression,  i.e., 

those  energy 

terms  due 

t°  Cj.  ^2. 

end  £q. 

°i=  ~T" 

|  iT  [  i  .  i 

♦  -j-  N 

1  4  2 

]  ddn 

(31) 

O 


where  K,  M^,  Hz  are  symmetric  18x18  arrays  given  in  Eqns 
(A8-A10)  in  the  Appendix. 

Each  term  of  the  K,  ,  and  expressions  of  Eqns 
(A8-A10)  represents  nine  terms  due  to  the  summation 
convention  on  i  and  j.  Each  of  these  nine  terms  results  in 
an  18x18  array  that  is  premultiplied  by  a  scalar  elasticity 
element,  see  also  Eqn  (22).  Eqns  (32)  show  one 
representative  term  from  K,  M ,  B„t,  respectively . 
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A11 

J 

0 

_ i 

j 

lxl 

18x1 

1x18 

A11 

.  °L1 

dT  o"l  ] 

(32) 

lxl 

18x1 

1x18  18x18 

»  :  A  [  H  d  dT  H  1 

2  11  L  °  1  0iJ 

lxl  18x18  18x1  1x18  18x18 

As  explained  by  [61],  the  desired  formalism  of  Eqns 

(25.26.28) ,  will  not  be  exhibited  using  the  terms  as  defined 

in  Eqns  (31 ,32 , A8-A10) .  That  is,  if  the  first  variation 
were  carried  through  on  Eqn  (31)  with  the  K,  B^,  and  Nz  as 
defined  in  Eqns  (A8-A10) ,  totally  new  terms  would  result. 
Similarly,  the  formation  of  the  linearized 

incremental/iterative  equations  would  result  in  still 
different  terms.  To  avoid  these  tedious  manipulations,  we 
seek  equivalent  energy  forms  for  K,  ,  and  Hz  such  that 
repeating  symmetric  terms  are  present  in  the  three 
expressions  for  n  ,  F(q),  and  K  ,  i.e.,  Eqns  (25,26.28). 

P  T 

Each  of  the  energy  terms  of  are  of  the  form  shown  on 
the  left  hand  side  (LHS)  of  Eqn  (33)  where,  refers  to 

the  elements  from  any  of  the  elasticity  arrays,  A_  to  j , 
p.r  =  0  to  7  and  i ,1=1,2, 6  as  before.  As  suggested  by  Bef 
[61],  if  each  of  these  terms  is  replaced  by  the  terms  of  the 
BHS  of  Eqn  (33) ,  the  desired  repeating  formalism  of  Eqns 

(25.26.28) ,  will  apply. 
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~~  fdTC.  .  L.  dT  H  d  d Ci  = 

2  J0  1J  P  i  r  j 


— i—  fdTC. .  [  L.  dT  H .  +  dT  L.  H.  ♦  H.  d  L*  Id 
6  |  i  j  Ip  i  r  j  pirj  pi  p  J  I 


dO 


(33) 

Using  the  underlined  term  from  of  Eqn  (A9)  with  the 
Eqn  (33)  substitution,  the  first  variation  is  as  shown  in 
Eqn  (34) .  It  is  stressed  that  this  term  is  representative 
of  all  the  terms  of  Eqn  (A9) . 


|  -i~  fdTB.  .  L.  dT  H.d  dO  1  = 

1  2  J„  "  °  1  *  J  J 

5  R  1: 


T  T 

L  d  H  +  d  L 
l  i  j  o 


i  .HJ  *  ,Bi  d  oLJ  ]d  dn  } 


— -  J  B  Jf(5d  L  d  H.d  ♦  d  L.6d  H  d  ♦  d  Ld  H  6d1 
eJuHoitJ  oi  ij  oiijj 

O  '■ 

+  K-Vi.-i-  *  dT<5dToLiiHid  +  dTdToLiiVd] 

+  (-TlWJ  d  ♦  dTiHi6d0Ljd  ♦  dTiHld0Lj  Sd  ]|  dO 


_l_  r 


B. 


J  "iJ 


O 


id’  [; 


3  L.d  H , 
o  i  t  j 


3,HldoL 


I] 


d  dfi 


(34) 


A  comparison  of  the  terms  in  the  second  and  final  lines  of 
Eqn  (34)  shows  that  the  variation  of  the  BHS  of  Eqn  (33) 
gives  a  repeating  form.  The  final  result  of  Eqn  (34)  is 
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expanded  into  unabridged  form  in  Eqn  (35) .  In  a  similar 
fashion,  ik  can  Shorn  that  the  derivation  of  the  linearized 
incremental/iterative  equations  (Eqn  (28))  also  results  in 
similar  repeating  forms. 


T-  J  {  B„  [w>i  *  ]“ 

4  BU  6dT  [3oLld\BJ  4  4  3,Vo^  ]< 


♦  B„„  «5d 
66 


t  r  t 

3  Ld 

1_  o  e  i 


H_  +  3d  L_  H_  ♦  3 
6  o  ot  6 


,v.1!  ]d } 


(35) 


Each  of  the  terms  of  the  X  and  arrays  are  of  the 
general  forms  of  the  LHS  of  Eqn  (36) .  In  Eqn  (36) ,  tmo 
forms  for  both  K  and  *2  are  shown.  With  the  once  again 
referring  to  elements  from  any  of  the  elasticity  arrays,  the 
LHS  are  those  energy  terms  from  Eqns  (3V.A8.A10',  tbs*.  «Io 
not  repeat  in  deriving  the  equilibrium  and  the  linearized 
incremental/iterative  equations.  If  each  is  replaced  by  the 
corresponding  RHS  of  Eqn  (36)  then  the  formalism  of  Eqns 
(25,26,28)  is  easily  shown  similarly  to  the  example  of 
Eqns  (33,34,35).  The  first  substitutions  shown  in  Eqn  (36) 
for  each  of  the  arrays  K  and  Xz  are  similar  to  those  of  Bef 
(61)  and  the  second  substitutions  are  unique  generalizations 
thereof . 


32 


Frow  K: 


+ Jr 


°ij  pLi  pLj  d  d°  (no 


on  p)  =  no  change 


-Hr 


2C.  L  L  d  dO  = 
ij  p  i  r  j 


-T-  \f  <=U  [  pN  pLj  ]<■  -O 


(p^r) 


Proe  H  : 
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■*  r  l  p  j  pi  Pj  2 


T  H.d  H . I 
P  1  p  J  J 


d  dO 


(p*r) 


Baaed  on  Eqna  (33,36),  Eqn  (31)  can  now  be  rewritten  as  in 
Eqn  (37) . 


d  dO 


r  t 

—  1  ~  1 

d* 

K  ♦  11  ♦  ~  h 

j 

*  4  2 

o 

J 

f  dT 

K  ♦  ♦  »z  1 

a 

L  3  6-1 
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where  the  elements  of  the  arrays  K,  ,  and  N,  are  found 
from  the  K.  ,  and  N2  of  Eqn  (A8-A10)  by  making  all  of  the 
aforementioned  substi tutions  of  Eqns  (33,36). 

The  transverse  shear  energy  terms  of  Eqns  (18,23,24) 
are  handled  in  a  totally  analogous  manner  only  much  more 
simply  since  transverse  shear  strains  are  linear  in 
displacement  and  hence  any  H  termis  are  zero  arrays  and 

J  IB 

therefore  not  needed.  Analogous  to  Eqn  (29) ,  the  transverse 
Shear  strains  of  Eqn  (15)  are  written  below  in  Eqn  (38) . 


s  T 

*  =  S  d 

m  o  m 


-  S  d 

2m  2  m 


(38) 


where  S  and  S  (m=4,5)  are  18x1  arrays  similar  to  the  L 
o  m  2  m  j  i 

of  Eqn  (29) . 

Due  to  the  assumed  linearity  of  these  strains,  there 

are  additional  energy  terms  from  the  transverse  shear 

strains  only  for  K  of  Eqn  (31) .  These  additional  terms  are 

found  from  Eqns  (23,24,38)  and  are  shown  below  in  Eqn  (39). 

The  equivalent  repeating  forms  that  are  added  to  K  of  Eqn 

(37)  are  found  via  Eqn  (36)  as  with  the  in-plane  energy 

terms  only  now  using  the  S  arrays  instead  of  the  L 

J  Bk  J  i 

arrays . 


where  m.n 


(A  S  ST  +  2D 
mn  o  m  o  n  mn 

4,5  as  before. 


S  sT 

o  ns  2  n 


F 

mn 


t 

S  S 
2  m  2  n 


)  d  dO 

(39) 


The  finite  element  is  defined  by 
interpolation  functions  used  to  approximate 


the  specific 
the  elements  of 
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this 


end  , 


the  displacement  gradient  vector,  d.  To 
approximate  the  continuum  displacements  by  shape  functions 
and  the  nodal  displacements  in  Eqn  (40) . 

u  =  Jf  q  (40) 

u  is  the  vector  of  continuum  displacements,  Jf  is  an 
array  of  interpolation  or  shape  functions,  q  are  the  nodal 
values  of  displacements  as  before.  Hext ,  the  displacement 


gradient  vector,  d 

,  is 

approximated 

using  Eqn 

(40)  as 

shown 

below  in  Eqn  (41) . 

d  =  h  q 

(41) 

where  X>  is  the 

array  of  shape 

functions 

and 

their 

derivatives.  Substituting  Eqn 

(41)  into 

the 

final 

expression  of  Eqn 

(37) 

including  the  transverse  shear 

terms , 

then  finally  gives 

Eqn 

(25)  for  the 

potential 

energy , 

where 

in  Eqn  (42)  we  can 

now 

identify. 

K  =  2)T  K  £  dO 

•  "i 

=  dt  a 

dfi.  »,  ■  1 

D  d~ 

O 

J 

o 

J 

1 

0 

(42) 

Vote  that  and  are  actually  functions  of  the  elements 
of  d  not  q  (see  Eqn  (36))  so  that  in  practice,  Eqn  (41)  is 
employed  prior  to  Eqn  (28)  in  the  solution  process.  The 
terms  of  Eqns  (25,26,28)  are  now  defined,  except  for  £> 
which  depends  on  the  specific  elemental  definition  and  will 
be  discussed  in  the  next  section. 

The  stiffness  arrays.  K,  ,  and  M2  are  formed  based  on 
literally  hundreds  of  18x1,  18x18,  and  1x18  symbolic  matrix 
multiplications,  see  Eqn  (A8-A10)  with  Eqns  (33,36) 
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substitutions.  These  multiplications  are  impractical  to 
carry  out  by  hand;  consequently.  the  symbolic  manipulator 
code,  MACSYMA,  was  employed.  Additionally,  MACSYMA  can 
convert  the  resulting  expressions  into  FORTRAN  statements. 
Since  these  arrays  are  independent  of  element  definition, 
subroutines  are  very  easily  implemented  that  calculate  the 
entries  of  K,  N^ ,  and  N2  where  elemental  information  is 
included  separately  through  £>. 

In  the  finite  element  solution  procedure,  the  stiffness 
of  each  element  is  calculated  individually  and  added  into  a 
global  stiffness  array.  In  this  way,  the  result  of  Eqns 
(25,42)  equivalently  represents  the  potential  energy  of  an 
individual  element,  where  the  total  n  is  found  by  summing 
the  energies  from  each  element. 

36  Degree  of  Freedom  Curved  Cylindrical  Shell  Element 

The  results  from  above  are  next  used  in  the  definition 
of  a  curved  cylindrical  shell  finite  element.  The  potential 
energy  expression  for  the  shell,  Eqn  (17)  with  cylindrical 
shell  special izations  for  the  strain  quantities,  Eqns 

.  o  1 

(A2-A6)  ,  requires  C  continuity  for  u,  v,  y .  and  a,*d  C 

continuity  for  w  [621.  The  simplest  approach  would  be  a 

curved  shell  element  derived  from  n  with  bilinear 

P 

interpolation  for  u,  v,  y^  ,  and  y^  and  cubic  Hermitian 
interpolation  for  w,  w, ^ ,  and  w , ^  resulting  in  a  4  noded  28 
degree  of  freedom  (dof)  curved  rectangular  element. 
However,  the  in-plane  continuum  displacement  dof,  u  and  v, 
are  coupled  to  the  transverse  displacement  dof  to  a  much 
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r  7*  ' 


7^ 


greater  degree  in  the  cylindrical  shell  compared  to  a  flat 
plate  due  to  the  nonzero  curvature  in  the  circumf erential 
direction  of  the  former.  Indeed,  in  a  linear  analysis,  u 
and  v  are  decoupled  from  the  remaining  dof  for  a  flat  plate. 
Consequently,  it  is  felt  that  linear  approximations  for 
these  variables,  though  quite  adequate  for  a  flat  plate,  are 
too  simple  for  the  shell,  especially  for  cases  where 
membrane  action  is  significant.  Therefore,  assuming 
quadratic  approximations  for  the  in-plane  dof  u  and  v,  and 
adding  four  nodes  located  at  the  midsides  of  the  element, 
results  in  an  eight  noded  36  dof  element,  see  Figure  4. 

Although  an  elesient  with  this  function  assumed  for  w  is 
nonconforming,  it  will  pass  the  patch  test  and  therefore 
converges  as  the  mesh  is  refined  (631.  However,  monotonic 
convergence  from  a  stiffer  solution  can  no  longer  be 
expected  as  in  a  conforming  elemental  formulation.  This 
incompatibility,  however,  is  not  nearly  as  significant  for 
the  present  formulation  as  it  is  for  simpler  ones,  e.g., 
Kircbhoff  plate  elements.  Recall  the  present  case  includes 
parabolic  transverse  shear  distribution  through  the 
thickness.  For  example,  the  shear  rotation,  ,  is  equal  to 
a  combination  of  w, ^  and  the  bending  rotation,  Wj .  see 
Figure  3.  For  the  limiting  thin  plate  and  shell  situation, 
we  would  expect  zero  shear  rotation  and  the  incompatibility, 
in  effect,  disappears  since  w, where  use  compatible 
approximations.  This  is  not  true  for  Kirchhoff  elements 
with  nonconforming  w  approximations.  It  should  be  pointed 
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development  [54]. 


/ 


FIGURE  4.  36  dot  Rectangular  Shell  Element, 

have  7  dof ,  Midside  Modes  have  2  dot. 


Corner  Modes 
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IV.  Lamlnatct.  Shell  In  Cvl  lndrleal  Bend  1  ng 


In  linear  analyses,  the  arrays  Hi  and  V2  of  Eqn  (26) 
are  not  used  and  the  equilibrium  equations  become  those 
given  in  Eqn  (43) .  Eqn  (43;  represents  a  system  of 
simultaneous  algebraic  equations  in  the  nodal  unknowns,  q, 
and  are  solved  by  gaussian  elimination.  All  elemental 
stiffnesses  are  integrated  exactly  usi.g  4x4  point  gaussian 
quadrature.  Upon  solving  Eqn  (43)  for  q,  the  elemental 
stresses  are  calculated  at  the  outermost  gauss  points,  i.e., 
those  nearest  the  corners  of  each  element.  The  stresses  are 
not  extrapolated  to  give  the  nodal  values. 

K  q  =  R  (43) 

The  present  theoretical  approach  and  finite  element 
casting  is  applied  to  a  laminated  circular  cylindrical  shell 
in  cylindrical  bending,  see  Figures  5  and  6  for  geometry  and 
material  properties.  The  shell  has  infinite  length  in  the 
x-direction;  is  simply  supported  at  the  two  ends,  s=0  and  L; 
and  is  subjected  to  a  sinusoidal  normal  pressure  loading. 
The  direction  of  the  pressure,  q(x,s),  will  induce 
compressive  stress  (c^)  on  the  bottom  of  the  shell  (C>0)  and 
tensile  stress  on  the  top  (C<0). 

The  infinite  shell  geometry  and  simple  supports  are 
represented  by  the  boundary  conditions  of  Eqn  (44). 
Additionally,  due  to  symmetry,  only  one  half  of  the  shell 
length  circumf erentially  is  modelled  and  the  halves  are 
matched  by  the  boundary  conditions  at  s=L/2.  The  length. 
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B,  shown  in  Figure  6  is  determined  such  that  all  elements 
have  an  aspect  ratio  of  one.  The  pressure  loading  is  given 
in  Eqn  (45)  . 


{ 


FIGURE  5.  Cylindrical  Shell  Geometry.  Radius  R=10*, 
Arclength  L- 10. 472',  thickness  h. 


FIGURE  6.  Planform  of  Shell  Showing  one  half 

Circumferential  Length,  L/2 ,  and  Discretization  Length,  B. 


along 

x  =  0  ,  B 

u 

=  W.J 

=  =  o 

(44) 

8  =  0 

w 

=  W.J 

=  0 

s  =  L/2 

V 

-  w’2 

=  =  0 

q  (x  ,  S  ,  0)  =  q 

sin 

( n s/L) 

(45) 

Resul ts 

are  generated 

f  or 

[90] 

unidirectional 

and 

(00/0/90]  cross  ply  laminates  and  compared  to  exact 
elasticity  solutions  given  by  Ren  (591.  Calculated 
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displacements  and  stresses  are  nona imena ional i zed  as  shown 
in  Eqn  (46).  The  transverse  displacement,  w,  and  the 
inplane  stresses,  o ^  and  ,  are  taken  from  the  center  of 
the  laminate,  i.e.,  at  S=L/2,  The  transverse  shear  stress, 
o-^ ,  and  the  circumferential  displacement,  ,  are  taken  from 
the  left  side  of  the  laminate,  i.e.,  at  s  =  0.  The  sign  of 
the  quantities  of  Eqn  (46)  is  opposite  that  given  by  Ben  in 
some  cases  since  the  respective  coordinate  systems  are 
different.  The  changes  are  made  so  that  the  physical  and 
nond i mens i ona 1 i zed  quantities  have  the  same  sign. 

Begular  meshes  of  5,  10,  and  20  elements  give  the 

nond i mens ional i zed  transverse  displacement,  w ,  for  S=100, 

[90]  laminate,  shown  in  Table  1.  Also  shown  is  the  exact 

value.  The  Table  shows  that  the  1x20  mesh  (x  by  s)  results 

converge.  Using  the  1x20  mesh,  further  displacement  and 

stress  results  are  generated  and  shown  in  Tables  2-9.  Two 

values  are  given  for  cr  in  Tables  5  and  9.  The  transverse 

4 

shear  stress  varies  very  rapidly  near  the  support  for  the 
thinner  shell  geometries  and  therefore  to  get  an  accurate 
value,  a  mesh  refinement  becomes  necessary.  The  value  in 
the  Tables  indicated  by  an  asterisk  is  the  stress  calculated 
based  on  a  1x20  mesh  that  has  smaller  elements  nearer  the 
support.  As  can  be  seen,  this  generally  gives  a  better 
value  of  o ^  for  the  thin  cases  but  is  unchanged  for  the 
thicker  laminates  where  the  shear  stress  has  a  smaller 
gradient  near  the  support. 
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S=B/h, 


(46) 


2=C/h 


Mesh 
(x  by  s) 

w 

Exact 

(Ren) 

1x5 

.0727 

1x10 

.0740 

1x20 

.0751 

.0755 

TABLE  1.  Convergence  of  Transverse  Displacement,  w,  for 
[00]  Laminate,  S=100. 


S 

W 

Exact 

(Ren) 

CST 

2 

.803 

.909 

.0764 

4 

.  278 

.312 

.0752 

10 

.  108 

.115 

.0749 

50 

.0762 

.0770 

.0748 

100 

.0751 

.0755 

.0748 

500 

.0746 

.0749 

.0748 

TABLE  2.  Transverse  Displacement,  «,  for  Various  Shell 
Thicknesses,  [90]  Laminate.  Classical  Solution  (CST)  taken 


""I" 


T7T 


7^ 


s 

°2 

Exsct 

(Ben) 

h/2 

-h/2 

h/2 

-h/2 

2 

-2.52 

1.336 

-2.455 

1.907 

4 

-1.219 

.950 

-1.331 

1.079 

10 

-  .  S36 

.773 

.890 

.  B  ,‘\> 

50 

-  .761 

.744 

-  .767 

.752 

100 

-  .759 

.742 

-  .758 

.751 

500 

-  .  777 

.718 

-  .752 

.750 

TABLE  3.  Circumferential  Stress,  ,  for  Various 

Thicknesses,  (901  Laminate,  for  (=h/2  (bottom)  and 
(top)  . 


s 

*1 

Exact 

(Ben) 

h/2 

-h/2 

h/2 

-h/2 

2 

- .0252 

.01336 

- .0245 

.0816 

4 

-.0122 

.00950 

-.0133 

.0260 

10 

-.00839 

.00773 

-.0089 

.0105 

50 

- .00761 

.00744 

-.0077 

.0076 

100 

-.00759 

.00740 

-.0076 

.0075 

500 

-.00777 

.00718 

-.0075 

.  0075 

TABLE  4.  Longitudinal  Stress,  ,  for  Various 

Thicknesses,  [90]  Laminate,  for  (,'  =  h/2  (bottom)  and 
(top)  . 


Shell 

C=-h/2 


Shell 

C=-h/2 
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s 

°4 

Exact 

(Bcr> 

2 

ft 

.489 

.489 

.555 

4 

ft 

.543 

.543 

.572 

10 

ft 

.558 

.559 

.570 

50 

ft 

.  528 
.557 

.  568 

100 

ft 

.449 

.543 

.  565 

500 

ft 

.  131 
.314 

.563 

TABLE  5.  Transverse  Shear  Stress,  a ,  for  Various  Shell 

Thicknesses,  [90]  Laminate,  for  £=0.  Values  with  an 
Asterisk  Result  from  Refined  1x20  Mesh. 


S 

W 

Exact 

(Ren) 

CST 

2 

1.141 

1.436 

.0799 

4 

.382 

.  457 

.0781 

10 

.  128 

.  144 

.0777 

50 

.0796 

.0808 

.0776 

100 

.0781 

.0787 

.0776 

500 

.0774 

.0773 

.0776 

TABLE  0.  Transverse  Displacement,  w,  for  Various  Shell 
Thicknesses,  [90/0/80]  Laminate.  Classical  Solution  (CST) 
taken  from  [59], 
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h/2 


-h/2 


Q 


Exact 

(Ren) 

h/2  -h/2 


2 

-3. 153 

1.732 

-3.467 

2.463 

4 

-1.406 

1.117 

-1.772 

1.307 

10 

-  .889 

.829 

-  .995 

.897 

50 

-  789 

.774 

.798 

.782 

100 

-  .787 

.770 

-  .786 

.781 

500 

-  .806 

.745 

-  .780 

.768 

TABLE  7.  Circumferential  Stress,  o-^ ,  for  Various 

Thicknesses,  [90/0/90]  Laminate,  for  [=h/2  (bottom) 
C=-h/2  (top) . 


S 

*1 

h/2 

-h/2 

Exact 

(Ren) 

h/2  -h/2 

2 

-.0316 

.0173 

-.0347 

.0871 

4 

- .0141 

.0117 

-.0177 

.0293 

10 

-.00889 

.00829 

-.0100 

.0115 

50 

-.00787 

.00774 

-.0080 

.0079 

100 

- .00787 

.00770 

-.0079 

.0078 

500 

-.00806 

.00745 

-.0078 

.0077 

TABLE  8.  Longitudinal  Stress,  ,  for  Various 

Thicknesses,  [90/0/90]  Laminate,  for  (-h/2  (bottom) 
O-h/2  (top). 


Shell 

and 


Shell 

and 
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s 

uT 

Exact 

(Ren) 

2 

.287 
*  .287 

.394 

4 

.326 
*  .326 

.476 

10 

.339 
*  .340 

.525 

50 

.  328 
«  .  340 

.  526 

100 

.  291 
*  .334 

.  523 

500 

.0853 
«  .217 

.  525 

TABLE  9.  Transverse  Shear  Stress,  o  ,  for  Various  Shell 

Thicknesses,  [90/0/90]  Laminate,  for  C=0.  Values  with  an 
Asterisk  result  from  Refined  1x20  Mesh. 

From  Table  2,  it  is  evident  that  by  allowing  transverse 
shear  deformation  in  the  unidirectional  laminate  the  shell 
response  is  much  more  flexible  in  the  thicker  laminates 
(compare  present  and  exact  to  the  CS'r  solution)  .  For 
S=R/h=10,  the  present  approach  underpredicts  the  transverse 
displacement,  w,  by  6X  yet  the  classical  solution  is  35% 
below  the  exact  result.  Similarly,  for  S=4,  the  present 
value  is  1 IX  below  the  exact  yet  the  CST  solution  is  76X 
less . 

For  thickness  ratios  of  S  =  4  and  lower,  we  cannot  expect 
close  agreement  between  the  present  2-D  approach  and  the 
exact.  Recall  that  the  theory  presented  in  Chapter  II 
clearly  uses  approximations  that  invalidate  the  approach  for 
the  very  thick  geometries.  Indeed,  for  the  S=  2  values  shown 
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in  the  Tables,  the  geometry  is  one  of  a  3-D  solid,  certainly 
not  a  thin  shell  (recall  earlier  as s limp t ions)  .  For  these 
cases  the  three  dimensional  effects  of  nonzero  transverse 
normal  stress,  o  and  a  varying  transverse  displacement 

u 

through  the  thickness,  w<'')  play  a  mr.jor  role  in  the  shell 
response.  Therefore,  the  entries  in  a1 '  Tables  for  S=2  and 
to  some  extent  S=4  as  well,  are  shown  for  completeness  only, 
although  the  majority  of  the  comparisons  to  the  exact  values 
are  still  quite  good. 

In  Table  2,  as  the  shell  becomes  very  thin  for  S=50  and 
above,  the  three  approaches  converge  to  the  same  result  as 
expected  since  the  transverse  shear  deformation  is  minimal 
and  classical  assumptions  are  valid.  See  also  Figure  7  for 
graphical  presentation  of  w  for  [90]  laminate. 

Similar  remarks  can  be  made  for  the  cross  ply  [90/0/90] 
laminate,  see  Table  6  and  Figure  8.  For  S=10,  the  present 
theory  gives  the  maximum  transverse  displacement,  w.  Ill 
below  the  exact;  the  classical  solution  is  46%  below  the 
exact.  Both  2-D  approaches,  i.e.,  present  and  CST,  are  in 
slightly  worse  agreement  compared  to  the  exact  than  in  the 
unidirectional  cases.  The  previous  arguments  for  the  [90] 
laminate  apply  to  the  cross  ply  laminate  in  addition  to  one 
other  comment.  The  assumed  displacements  give  continuous 
strains  with  respect  to  the  transverse  coordinate,  C- 
However,  discontinuities  in  stress  result  due  to  different 
constitutive  relations  in  the  cross  ply  laminate  for  each 
ply.  Although  elasticity  theory  also  predicts  these 
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discontinuities  only  for  the  inplane  stresses,  the 
transverse  shear  stress,  ,  must  be  continuous  so  that 
equilibrium  is  satisfied.  The  present  2-D  approach 
satisfies  equilibrium,  generally,  only  at  the  midplane  of 
the  shell.  The  consequence  of  this  becomes  evident 
comparing  Tables  2  and  6,  the  transverse  displacement 
response  for  [90]  and  [90/0/90]  laminates.  Further 
consequences  are  found  by  examining  the  stress  results. 

Tables  3  and  7  show  the  results  for  the  circumf erential 
stress,  cr^  for  both  laminate  types.  As  shown,  this  stress 
is  calculated  very  well  for  both  the  unidirectional  and 
cross  ply  laminates  even  for  the  thick  cases.  Recall  that 
the  stresses  shown  for  the  present  case  are  found  at  the 
outermost  gauss  points  within  each  element.  These  points  do 
not  correspond  exactly  to  the  center  of  laminate  as  the 
Table  indicates.  Increased  accuracy  in  o ^  could  be  attained 
by  a  mesh  refinement  near  the  center  of  the  shell  thereby 
smving  the  gauss  points  closer  to  the  actual  shell  center  or 
by  extrapolating  the  stresses  to  give  nodal  values. 

Improving  the  numerical  aspects  of  the  solution  as  just 
described  will  increase  the  accuracy  somewhat  for  the 
thicker  shell  geometries  as  well.  However,  the  largest  loss 
of  accuracy  for  these  cases  is  due  to  the  2-D  assumption  of 
the  present  approach,  i.e.,  neglect  of  o',  in  the  solution. 
Notice  in  Tables  3  and  7  that  the  present  results  are  the 
most  inaccurate  at  the  top  of  the  shell,  i.e.,  where  the 
load  is  applied,  and  therefore,  where  the  effects  of  ct 

3 


are 


generally  the  moat  significant.  The  circumferential  stress, 
o  ,  is  also  clotted  as  a  function  of  thickness  in  Figures 
9-12. 

Tables  4  and  8  give  the  or.gi  iudinal  stresses,  ,  in 

the  shell.  The  comparisons  for  bot-  laminate  types  to  the 
exact  are  fairly  good  except  for  thicknesses  of  S=4  and 
thicker.  In  the  very  thick  cases,  the  longitudinal  “tress 
is  significantly  underestimated  on  the  top  surface  of  the 
shell,  i.e.,  at  C=“h/2.  Apparently,  the  influence  of  the 
three  dimensionality,  i.e.,  cr  on  the  top  of  the  shell  where 

sJ 

loading  acts,  has  a  great  effect  for  this  stress.  The 
influence  is  much  less  at  the  bottom  of  the  shell  where 
is  identically  zero,  see  C *59 1  . 

The  transverse  shear  stress,  a  is  calculated  at  the 
midplane  of  the  unidirectional  laminate  very  well  for  values 
as  small  as  S=4 ,  see  Table  5.  The  stress  for  S=500  as 
predicted  in  the  present  approach  is  inaccurate  because  its 
gradient  is  so  high  for  this  thickness  near  the  support.  A 
more  accurate  value  should  be  attainable  if  a  further 
refined  mesh  were  analyzed  and  theeref ore ,  the  gauss  point 
where  is  calculated  would  be  closer  to  the  end  support. 
Figures  13-15  show  the  transverse  shear  stress  as  a  function 
of  thickness  for  the  unidirectional  case.  Notice  that  as 
the  laminate  becomes  thicker,  the  exact  transverse  shear 
stress  becomes  more  asymmetrical ,  whereas  the  present 
approach  always  gives  a  symmetric  response  about  the 
midplane.  Becall  that  the  cubic  term  in  the  stress 
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expression  for  the  present  approach  which  gives  this 
asymmetry  was  neglected  as  small,  see  text  following  Eqn 
(16).  The  Figures  show  that  this  is  a  good  assumption  since 
the  asymmetry  is  significant  only  for  the  thick  S  =  4  case. 

Similarly,  the  transverse  stress,  a  ,  is  given  for  the 
cross  ply  laminate  in  Table  9  and  in  Figures  16-18.  As 
argued  previously,  the  discontinuities  in  as  a  function 
of  transverse  coordinate,  £,  violate  equilibrium  away  from 
the  shell  midplane  and  thus  explains  the  disagreement  shown. 
Hotice  that  the  results  from  the  present  theory  at  the  ply 
interfaces  are  approximately  equidistant  from  the  exact 
value  at  the  interfaces. 

Finally,  Figures  19  and  20  show  the  circumferential 
displacement,  V,  taken  at  the  left  support  for  both  laminate 
types.  The  present  approach,  with  an  assumed  cubic  function 
through  the  thickness,  Slightly  underestimates  the  exact 
values  and  but  captures  the  distribution  very  well.  At  the 
midplane  where  is  typically  maximum,  the  shear  deformable 
theories  give  a  slightly  steeper  slope  compared  to  the  CST 
result  and  gradually  approach  the  same  slope  at  the  outer  C 
coordinates  of  the  shell  where  the  shear  deformation  is 
zero . 
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w 


S  =  R/h 


FIGURE  8.  Transverse  Displacement,  w,  versus  Shell 


Thickness,  S,  for  [90/0/90]  Laminate. 


FIGURE  8.  Circumferential  Stress,  ,  Through 

Thickness  for  S  =  4,  [90]  Laminate. 
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•  Present 


0.50  -> 


FIGURE  10.  Circumferential  Stress, 
Thickness  for  S=10,  [90]  Laminate. 


Through 


FIGURE  11.  Circumferential  Stress.  .  Through 

Thickness  for  S=4 ,  [00/0/00 ]  Laminate. 


FIGURE  12.  Circumferential  Stress,  Through 

Thickness  for  S=10,  [90/0/901  Laminate. 


FIGURE  13.  Transverse  Shear  Stress 
Thickness  for  S-4,  [90]  Laminate. 


FIGURE  14.  Transverse  Shear  Stress,  o  ,  Through 


the 


Thickness  for  S~10,  [SO]  Laminate. 


0.20 


•  Present 


FIGURE  IS.  Transverse  Shear  Stress,  ,  Through  the 

Thickness  for  3=100,  [90]  Laminate. 


0.50 


•  Present 


FJQUBE  17.  Transverse  Shear  Stress,  ,  Through  the 
Thickness  for  S  = 1 0 ,  [90/0/00]  Laminate. 


Thickness  for  S=100,  [90/0/90]  Laminate. 


FIGURE  20  Circumferential  Displacement,  V,  Through  the 
Thickness  for  S=10,  [90/0/90]  Laminate. 
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V .  Conclua iona 


A  higher  order  8hell  finite  element  theory  haa  been 
developed  that  Includes  a  through  the  thickness  parabolic 
shear  stress  distribution.  The  element  has  36  degrees  of 
freedom  and  avoids  the  locking  phenomenon  usually  associated 
with  shell  elements.  Laminated  compo»if e  construction  can 
be  considered.  Exact  Green  strain  rel^uons  are  included  so 
that  large  displacement,  simplified  large  rotation  (up  to 
the  assumption  of  linear  strain  displacement  relations  for 
tiie  transverse  strains)  degrees  of  freedom  can  be 
investigated.  The  approach  is  applied  to  a  laminated 
cylindrical  shell  in  cylindrical  bending  and  linear  results 
are  compared  to  the  exact  elasticity  solution. 
Displacements  and  stresses  are  predicted  very  well  for  the 
unidirectional  laminate  even  for  very  thick  geometries. 
Stresses  are  predicted  less  accurately  for  the  cross  ply 
laminate  due  to  violation  of  equilibrium  away  from  the 
midplane  of  the  shell,  a  consequence  of  the  2-D  approach. 
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